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Nonlinear corrector for RANS equations 


Loic Frazza* Adrien Loseille^ Frederic Alauzetd Alain Dervieux § 
Sorbonne Universites, UPMC, 4 place Jussieu 75252 Paris cedex 05, France 
GammaS Team, Inria Saclay, 91120 Palaiseau, France 


The scope of this paper is to present a nonlinear error estimation and correction for 
Navier-Stokes and RANS equations. This correction is obtained by deducing a source term 
from the evaluation of the residual of the solution interpolated on the h/2 mesh. To avoid 
the generation of the h/2 mesh (which is prohibitive for realistic applications), the residual 
at each vertex is computed by local refinement only in the neighborhood of the considered 
vertex. It successfully improves solution predictions and yields a sharp estimate of the 
numerical error. 


I. Introduction 

The numerical simulation of engineering problems involves a discrete solution which is alleged to converge 
toward the continuous solution of the problem as the mesh-size decreases. As the exact analytical solution 
of the problem is sought, the difference between the numerical and the analytical solution can be seen as a 
noise. In an engineering context, this noise can have disastrous consequences on the numerical prediction 
which can, in turn, lead to actual accidents or misconceptions. 

It is thus essential to reduce and estimate this error. This is usually done by leading a grid convergence 
study, where the error is estimated with the difference between two solutions computed on two successives 
grids. However, this approach requires an additional finer grid (thus expensive) to appropriately estimate the 
error. Moreover, it requires the expertise of an engineer. In a mesh adaptation context multiple successive 
grids are automatically generated following a sensor which estimates the error related to the local mesh size. 
It does not requires any mesh prescription, but it needs a relatively accurate estimation of the numerical 
error to be efficient and to stop when the mesh and the solution are converged. 

Recently, different adjoint based mesh adaptation strategies 1-6 have been developed in order to optimize 
the mesh with respect to a specific functional output. These techniques uses the adjoint solution to estimate 
the error committed on the output and minimize it with respect to the mesh size. However, these approaches 
are based on a linearization of the equations and can thus be inefficient for strongly nonlinear problems like 
RANS equations. This is why we proposed a nonlinear equivalent of these approaches in order to efficiently 
compute a correction in the global solution and in the output functional (lift, drag, heat transfer, sound, ...). 

*PhD student, Sorbonne Universites, UPMC Univ Paris 06, Gamma3 Team 
t Researcher, Gamma3 Team 

Researcher, Gamma3 Team 
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II. Adjoint based error estimation and correction 


A. Problem setup and Numerical discretization 

1. Navier-Stokes RANS equations 

The compressible Navier-Stokes equations for mass, momentum and energy conservation read: 


g(pu) 
dt 
d(pE) 

dt 


^ + V • (pu) = 0, 

+ V • (pu ® u) + Vp = V-T, 

+ V-((p£7 + p) u) = V-(T-u) + V-(AVT), 


(1) 


where p denotes the density (kg/m 3 ), u the velocity (m/s), E the total energy per mass (m 2 .s -2 ), p the 
pressure (N/m 2 ), T the temperature ( K ), p the laminar dynamic viscosity (kg/(m.s)) and A the laminar 
conductivity. T the laminar stress tensor: 


T=(p + p t ) 


(V <g> u + 4 V <g> u) 


-V.ul 


where (in 3D) u = (u, v, w) and 


V.ul 



U X +Vy +W Z 


where u x = u y = u z = (idem for v and w). 

The variation of nondimensionalized laminar dynamic viscosity and conductivity coefficients p and A as 
a function of the dimensional temperature T is defined by Sutherland’s law: 


A 



{ Too + Su A 

V T + Su ) 


and 



/ Too + SuA 

U+Su ) ’ 


where Su = 110 is the Sutherland constant and the index oo denotes reference quantities. The relation 
linking p and A is expressed from the Prandtl laminar number: 

Pr = with Pr = 0.72 for (dry) air, 

A 

where C p is the specific heat at constant pressure. 


2. Turbulence modeling 

According to the standard approach to turbulence modeling based upon the Boussinesq hypothesis, the 
turbulence is modeled with an eddy viscosity pt, which is added to the laminar (or dynamic) viscosity p. 
The dynamic viscosity is usually taken to be a function of the temperature, whereas pt is obtained using 
a turbulence model. Here we choose the Spalart-Allmaras one equation turbulence model 7 given by the 
following equation: 


dv 

dt 


+ u • Vj> = cm [1 - fn]Sv - 



Cbl , 
o Jt‘2 


p i + ~ [ v ' i( v + t')Vz>) + c 62 ||Vf>|| 2 ] + /uAu 2 , (2) 


where v is the turbulent kinematic viscosity and all the constants are defined below. In the standard model 
the trip term is being left out, i.e., ft i = 0. Moreover, some implementations also ignore the fa term as it 
is argued that if the trip is not included, then fa is not necessary. 8 In Wolf, this simplified version has been 
considered and we prefer to write it under the following form, which is more appropriate for its discretization 
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with the finite element/finite volume method. Indeed, Equation (2) can be decomposed into the following 
terms: 

~K7~ + u ■ Vpd = c b iSpv -c w if w p(-') +-V-((i/ + i>)Vi>) + ^^||Vd|| 2 . 
at '—v—' '—v—' \a) a a 

convection production ^ ~~ v v 

destruction dissipation diffusion 


Notice that this is not a conservative model. If a conservative form of the Spalart-Allmaras is foreseen, we 
have to consider the variation proposed by Catris and Aupoix. 9 The turbulent eddy viscosity is computed 
from: 

lit = pvfvi, 


where 


fv 1 


r 

x 3 + c 


vl 


and 


v 



Additional definitions are given by the following equations: 


with v = — . 

P 


fv2 — 1 ~ 


X 

1 + xfv 1 


and 


S = fl + 


v 

n 2 d 2 


fv 2 


where = ||V x u|| . 


d is the distance to nearest wall which is computed for each vertex at the beginning of the simulation. The 
set of closure constants for the model is given by 


17 = 


2 

3 ’ 


Cwl 


Q)1 

K 


1 + 0,2 

<J 


Chi = 0.1355, 0,2 = 0.622, k = 0.41, 
Ot;2 — 0.3 , Cw3 2 , OjI — T.l . 


Finally, the function f w is computed as: 


/ 1 + c 6 , \ 1/6 

fw=9[ 6 , X ) with 9 = r + c w 2 (r 6 - r) 
\y ' ^w3J 


and 


r = min 


v 

Su 2 d 2 



5. Vector form of the RANS system 

We write the RANS system in the following (more compact) vector form: 

W t 4 F\{W) X 4 F 2 (W) y 4 F 3 (W) z = Si(W% + S 2 (W0*/ + 5 3 (W) z + Q(W), 

where 5*(W) 0 = - (■i = 1,2,3, a = x,y,z) (idem for F). W is the nondimensionalized conservative 


variables vector: 


da 


W = (p, pu, pv, pw, pE, pv) . 

F(W) = (Fi(W), F 2 (W), F 3 (W)) are the convective (Euler) flux functions: 



F\(W) 

= 

(pu, 

pu 2 d 

- p, puv, puw, u(pE + p), puv) T , 


F 2 (W) 

= 

(pv, 

puv , 

pv 2 + 

p , puw, v(pE + p), pvi>) T , 


F 3 {W) 

= 

(pw, 

puw , 

pvw , 

pw 2 4 

■p, w(pE + p), pwv) T 

),S 2 (W),S 3 (W)) 

are 

the laminar viscous fluxes: 

Si(W) 

= (°’ 

f~xx , 

f~xy, 


U T 

^ ' XX 

4 vFxy 

+ w71 z 4 AT/, — (v 4 

(7 

S 2 {W) 

= (°’ 

f~xyi 

fyy, 

fyz, 

U~Txy 

4 uT/y 

4 ivFyz 4 AT/, — (y 4 

S 3 (W) 

= (°’ 

T xz , 

fyz, 

Tzz, 

uTxz 

4 u7^ z 

4 wl~ zz 4 AT/, — (^ 4 i 

<7 


( 3 ) 


(4) 
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where 7 1j are the components of laminar stress tensor T = (p + pt) [(V (g> u + 4 V <8> u) — |V. ul] : 
2 

% x = (p + Ht) g (2Ms -v y -w z ), T xy = (p + p t ) ( u y +v x ), T xz = (p + p t ) ( u z + w x ), 


Q(W) are the source terms, i.e. the diffusion, production and destruction terms from the Spalart- 
Allmaras turbulence model: 

Q(W) = ( 0, 0, 0, 0, 0,^\\Vu\\ 2 + pc bl Su + c wl f wP (^j ) T . (5) 

Note that Q = 0 in the case of the laminar Navier-Stokes equations, unless additional source terms are added 
(to take into account gravity, for instance). 


4- Spatial discretization 

The spatial discretization of the fluid equations (1) and (2) is based on a vertex-centered finite element/finite 
volume formulation on unstructured meshes. The equations are integrated on each finite volume cell C) (using 
the Green formula): 

\C i \^+F i = S i + Q i , ( 6 ) 

at 

where IT) is the mean value of the solution W on cell C), F, : , S,; and Q, are respectively the numerical 
convective, viscous and source flux terms: 

F< = f F(Wi ) • iijdy, S, = f S(Wi) ■ n, d 7 , Q, = f Q(W,) dH, 

JdCi JdCi JCi 

where rq is the outer normal to the finite volume cell surface dCi , and F, S and Q are respectively the 
convective, viscous and source terms flux functions as defined previously in Relations (3), (4) and (5). The 
numerical scheme combines a HLLC approximate Riemann solver 10 to compute the convective fluxes and the 
Galerkin centered method to evaluate the viscous terms. Second order space accuracy is achieved through 
a piecewise linear extrapolation based on the Monotonic Upwind Scheme for Conservation Law (MUSCL) 
procedure 11 which uses a particular edge-based formulation with upwind elements. 


HLLC approximate Riemann SOLVER. The idea of the HLLC flow solver is to consider locally a simplified 
Riemann problem with two intermediate states depending on the local left and right states. The simplified 
solution to the Riemann problem consists of a contact wave with a velocity Sm and two acoustic waves, 
which may be either shocks or expansion fans. The acoustic waves have the smallest and the largest velocities 
(Si and Sj , respectively) of all the waves present in the exact solution. If Si > 0 then the flow is supersonic 
from left to right and the upwind flux is simply defined from F(Wi) where H) is the state to the left of 
the discontinuity. Similarly, if Sj < 0 then the flow is supersonic from right to left and the flux is defined 
from F(Wj) where Wj is the state to the right of the discontinuity. In the more difficult subsonic case when 
Si < 0 < Sj we have to calculate F(W*) or F(W*). Consequently, the HLLC flux is given by: 


foHLLC 


!U,.IT,.ti 0 S 


F(Wi) ■ n ij 

if 

Si > 0 

F(W*) ■ n ^ 

if 

Si < 0 < Sm 

F(Wj) ■ n ^ 

if 

Sm < 0 < S) 

F(Wj) ■ n ^ 

if 

Sj < 0 


W* and Wj are evaluated as follows. We denote by 77 = u ■ n. Assuming that 77 * = rf* = rj* = Sm, the 
following evaluations are proposed 10 (the subscripts and j are omitted for clarity): 


W* 


1 

S-S M 


p(S-v) \ 

pu (S - 77 ) + (p* - p)n where p* = p (S - r])(S M - 77 ) + p. 
pE (S — rj)+ p*S M — PV J 
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A key feature of this solver is in the definition of the three waves velocity. For the contact wave we consider: 

c = PjVj( s J - Vj) ~ PWiiSi - Vi)+Pi ~ Pj 
M Pj(Sj — Vj) ~ Pii^I — r]i) 

and the acoustic wave speeds based on the Roe average: 

Si = — Cj, 77 — c) and Sj = max(r]j + Cj,fj + c) . 

With such waves velocities, the approximate HLLC Riemann solver has the following properties. It automat¬ 
ically (i) satisfies the entropy inequality, (ii) resolves isolated contacts exactly, (iii) resolves isolated shocks 
exactly, and (iv) preserves positivity. 

2ND-ORDER ACCURATE VERSION. The MUSCL type reconstruction method has been designed to increase 
the order of accuracy of the scheme. 11 The idea is to use extrapolated values Wjj and W Jt instead of W, 

and Wj at the interface dCij to evaluate the flux with the approximate Riemann solver. Note that, in the 

implementation, the primitive variables (p, u, p ) are extrapolated to guarantee the positivity of the density 
and the pressure, then the conservative variables are reconstructed from these values. Thus, the gradients of 
the primitive variables are evaluated. However, in the following, we still denote by W the primitive variables 
vector. The numerical flux becomes: 

*ij = *? j LLC (W ij ,W ji) n ij ), 
where Wij and Wji are linearly extrapolated as: 

Wn = W Z + ^(\7W) Z -P^ and W jt = Wj+ ^(VW) j -P~Pi. 

In contrast to the original MUSCL approach, the approximate ’’slopes” (V1U)„ and (VlU)j,; are defined for 
any edge and obtained using a combination of centered, upwind and nodal gradients. 

The centered gradient, which is related to edge I\Pj , is implicitly defined along edge P t P 3 by the relation: 

(VWOg ■P^j = Wj-W i and (VW)g ■P~f i = W i -W j . 

Upwind and downwind gradients, which are also related to edge PiPj , are computed according to the 
definition of upwind and downwind tetrahedra of edge PiPj ■ These tetrahedra are respectively denoted K t j 
and I\ji . Kij (resp. Kj,) is the unique tetrahedron of the ball of Pi (resp. Pj) the opposite face of which is 
crossed by the line defined by the edge PiPj , see Figure 1. Upwind and downwind gradients are then defined 



Figure 1. Downwind Kij and upwind Kji tetrahedra associated with edge PiPj . 
for vertices Pi and Pj as: 

(VW0g = (VW)k, and (VW)g = (VW)\ Kji . 

where {VW)\k = J2paK WpV(I>p\k is the Pi-Galerkin gradient on tetrahedron K. Parametrized gradients 
are built by introducing the /3-scheme: 

(VW)* • PiPj = (1 — /3)(VW)^ • Pp^j + /3 (VW)Yj ■ pPj 

(yw)j-p^i = (i-p)(vw)fj-p^i + p(vw)P-p^i, 
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where /3 £ [0,1] is a parameter controlling the amount of upwinding. For instance, the scheme is centered 
for /3 = 0 and fully upwind for (3=1. 


Fourth-order numerical DISSIPATION: V4-SCHEME. The most accurate /3-scheme is obtained for 
(3 = 1/3. Indeed, it can be demonstrated that this scheme is third-order for the two-dimensional linear 
advection on structured triangular meshes. On unstructured meshes, a second-order scheme with a fourth- 
order numerical dissipation is obtained. These high-order gradients are given by: 

(VW)Y 4 -P^ = ? (VFF)g • i (VWOg • iy? 

? (VfF)g • iyi* + i (VW)% ■ p^i . 


The MUSCL schemes are not monotone and can be a source of spurious oscillations especially in the 
vicinity of discontinuities . 12 These oscillations can affect the accuracy of the final solution or simply end 
the computation because (for instance) of negative pressures. A widely used technique for addressing this 
issue is to guarantee the TVD property in ID 13 or the LED property in 2D/3D of the scheme, which ensures 
that the extrapolated values W tJ and IF/,; are not invalid. To guarantee the TVD or the LED properties, 
limiting functions are coupled with the previous high-order gradient evaluations. The gradient is substituted 
by a limited gradient denoted (VPF)^ m . The choice of the limiting function is crucial as it directly affects 
the convergence of the simulation. In this work, we use the Piperno limiters 14,15 which is expressed in a 
factorized form, 

(yW) Lim = VIT C $ PI ( VW C ^ 


\\7W V4 J ’ 


with 


i’pi(R) = (g + g-R) 


■li_«2_ 

J R 2 U R 

J_ _ Qi _ 
R 3 °R 


-19 


i+ 4s +i), s - i)3 


if R < 1 


if R > 1 


where R = 


\/W c 

vw V4 ' 


B. Adjoint based correction 

1. Linear correction: 

Despite some differences, the usual adjoint based approaches relies on the following implicit correction: 1-6 
given an analytical problem defined as 

'P(IV) = 0 , 

a numerical discretization of the problem is determined (using here finite volume approach) 

*h{W h ) = 0. 

From the numerical solution an output functional jh(Wh) is computed. Assuming the numerical solution 
Wh is close to the analytical solution W, the functional j can be linearize as 

j h (U h W) « j h {W h ) + J • (Iw - W h ), 

where II hW stands for the projection of the exact solution on the mesh and J = is the gradient of j 
with respect to the solution. The numerical equation can also be linearized as in the Newton’s method as 

*fc(IW) « * h (W h ) + A ■ {W h W - W h ), 

where A = ^ 7 - is the jacobian of Put together we obtain 

jhQihW) w j h (w h ) + w* ■ 9 h (n h w), 

where W* = J ■ A -1 = A~ T ■ J is the adjoint of the problem. This formulation is convenient for mesh 
adaptation as W* can be precomputed separately and 4/(11^ W) can be estimated with local interpolation 
errors. However, we can see that intrinsically this is equivalent to compute a corrected solution as 

W = W h + A~ 1 ^ h {U h W) 
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and then compute the functional j 

jh(W)*j h (W h ) + J-(W-W h ), 
so that the correction on itself is linear. 

Similarly, error transport equations are based on the linearization of the given problem: 16-19 a linear 
correction problem is deduced from the initial problem and thus provide a linear correction of the solution. 
The linearization of the problem is thus a determining factor of the quality of the corrected solution. 20 

2. Nonlinear correction: 

This is why we propose a nonlinear approach for finite volume solvers based on the fact that 

^ h (IL h W) =e h ^0 


and 

*(W h ) =e^0. 

In the finite volume approach, the residual is defined as the integral of the continuous residual of a recon¬ 
structed solution over a cell or equivalently for inviscid flow as the integral of the fluxes through the boundary 
of the cell. In that sense, for a given mesh an infinite family of continuous fields W verifies W) = 0, the 
reconstructed continuous field based on the discrete solution is one of them. 

However, if we consider finer meshes, in particular the n th subdivisions of the initial mesh, the only 
held verifying 'l' h/niW ) = 0 for any submesh is the solution of the continuous problem: W = W. Thus, 
the residual of the discrete solution interpolated on a finer mesh 'l' h/ 2 {Ih/ 2 Wh) = £h /2 gives us a usefull 
indication of the quality of the solution. 

We propose the following approach, inspired by multigrid methods, to compute a corrected solution Wh 
from the residual ^h/ 2 {Ih/ 2 ^h)- The mesh is initially divided by two (or a power of two) and the solution 
Wh is linearly interpolated on it. The residual is then computed on the finer mesh and accumulated back 
on the initial mesh. 

The linear accumulation process can be seen as an inverse operation of the linear interpolation: The 
interpolated value of a vertex of a finer mesh is a linear combination of the values of some vertices of the 
coarse mesh. In the accumulation process, the residual of this vertex of the fine mesh will be added to the 
contribution of the vertices of the coarse mesh weighted by the same coefficients. This forms a source term 

Sh = Ih/2->h (*/,/2 (h^h/2(W h ))) , 

which is then introduced in the equation to compute the correction as 

*h(W h ) = S h . 

The corrected solution Wh is computed performing a few iterations with the same solver, with the source 
term Sh, which nonlinearly propagates the errors. This approach automatically encompasses all features of 
the flow solver (order, flux reconstruction, limiters, ... ) both in the flux computation and the resolution. 
All of the operations involved here are local, so that we don’t actually need to generate the actual h/2 mesh, 
which would be far too expensive. Instead, we extract for each point its vicinity mesh and divide it by two 
to compute and accumulate the source terms fluxes. This approach yields no memory over cost and can be 
done efficiently in parallel. Finer subdivision (h/ 4, h/8, ...) of the initial mesh can also be used to improve 
predictions. The overall procedure is thus: 

• Compute the solution on the initial mesh: W/, 

• Linearly interpolate the solution on the h/2 mesh: lI/ l / 2 W / / i 

• Compute the h/2 residual: ’P h/ 2 {Ih/ 2 ^h) = £h /2 

• Accumulate the h/2 residual on the h mesh: Sh = A'^ 2 (E h / 2 ) 

• Perform a few iterations of the flow solver with Sh as source term: IF h(Wh ) = Sh- 
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Numerical Implementation: Dividing the element size by two multiplies the number of elements by 
4 in 2D and 8 in 3D, it is thus particularly expensive to fully generate the submeshes for large practical 
cases. As the residual computation and accumulation is a local operation, it can be done at a local level, by 
generating virtually the submeshes of the ball around each vertex. That way no additional memory cost is 
implied and the process can straightforwardly done in parallel. Similarly, multiple computations of the same 
element can be avoided with a careful ordering detailed here after. 

Figure 2 shows an example of 2D mesh and it’s subdivision by two. Starting with the solution on the 
coarse mesh (circle dots and dark lines), the solution is interpolated on the divided mesh (grey lines and 
square dots). The residual is computed on the fine mesh and accumulated on the coarse mesh (weighted by 
the barycentric coordinates). 

A careful analysis of the numerical discretization presented in Section II.A shows that the computation 
of the residual of a given vertex requires the values of the vertices of its third order ball (neighbours of the 
neighbours of its neighbours). Similarly, to compute the accumulated source term of a given vertex, we only 
need the vertices of the fine mesh encompassed by the first oder ball of the given vertex in the coarse mesh. 
This is why, we can generate separately the subdivided mesh of the second order ball for each vertex of the 
coarse mesh to compute its source term. 

However, it can be reduced to the subdivision of the first order ball if we consider carefully the different 
terms involved in the computation of the residual. In our vertex centered finite volume solver Wolf, the 
advection fluxes are computed and summed for each edge, while the viscous terms are computed for each 
triangle. We can thus compute them separately in any order as explained below. 

We can see in Figure 2 that the residual of the fine node C contributes through the accumulation process 
to the source term of coarse node A and B by a factor 1/2. The residual of node D contributes to the 
source term of node A but not of node B (barycentric coordinate 0), however, its solution is required to 
compute the residual of node C and thus contributes indirectly to the source term of B. Finally, node E 
directly contributes to the source term of A, but it belongs to an upwind element of the edge DC, so that it 
contributes to the residual of C and the source term of B. 

If we want to compute the source term of node B, the residual of node D does not need to be computed, 
and the residual of node C which contributes to both source terms of A and B can be computed only once. 
To do so, we can notice that the advection term in the residual of C is the sum of the fluxes across each 
edge connecting C to one of its neighbours, i.e. here the edges CA, CD, CB, ... As the residual of C is 
added to the source terms of nodes A and B in the accumulation, we can compute and add separately the 
contributions of edges CA and BC. 

To this end, we generate the local sub-mesh as shown in Figure 3. The grey part of the sub-mesh is 
not needed and not generated at all, the black part of the sub-mesh is generated for the interpolation of 
the nodal values and the fluxes are computed only on the red edges and triangles encompassed. That way, 
only a half of the residual of node C is computed and accumulated on both nodes A and B and the other 
half of the residual is computed in the vicinity of A. Additionally, the central sub-triangle contribution is 
computed in the vicinity of the node with the smallest ID in the coarse triangle. Each contribution of each 
fine edge and triangle are thus computed only once which reduces computational costs. It also ensures that 
each edge has an upwind and downwind triangles for second order extrapolations in the HLLC flux solver 
(see Section II.A.4). 


III. Numerical Results 


In the sequel, given the analytical solution of a problem W and the numerical solution Wh, we split the 
approximation error W — Wh into the interpolation error and the implicit error as 

W-W h = IF — 11/, H) + n fe W - Wh - 

Interpolation Error Implicit Error 

The interpolation error is entirely related to the discretization of the solution and can only be improved by 
a refinement of the mesh. Here, we are interested in the implicit error as we want to improve the discrete 
value of our solution. We can thus compute analytically the exact implicit error: 

£eZ = ^£\\n h w-w h \\*dn, 
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Figure 2. Mesh Subdivision 


Figure 3. Mesh used for the computation of the 
source term for node B 


the estimated implicit error: 

^ £ \\w h -w h \\*dn, 

as they can be expressed element-wise as polynomial expressions. Similarly, we define a function j computed 
from u , and we define the exact functional error: 

C = l)W)-i(ni.n 

the estimated functional error: 

££“" = \j(W h ) - j(W h )\, 

and the corrected functional error: 

& = \j(w h ) - j{n h w)\. 

Note that the interpolation error is not taken into account here as, we compare the numerical functional 
j(Wh) to j{JlhW) which is an approximation of the exact j(W). However this interpolation error can be 
taken into account in another means. In the sequel, we will use a finer solution as estimation of the exact 
solution. 

A. 2D Laminar NACA0012 

We first validate the present corrector in an adaptive context on a NACA0012 geometry, at a Reynolds 
number of 500 with a purely laminar viscous flow. The Mach number is fixed to M = 0.1 and the angle 
of attack to a = 3°. Adapted meshes are generated iteratively with feflo.a 21,22 using a goal oriented 
error estimate on the drag coefficient. 6 An example of the resulting meshes and solutions is shown in 
Figures 4 and 5. 

Figure 8 shows the convergence of the drag computed for both the initial (red) and corrected (blue) 
solution. We can see that corrector moderately improves the solution and the computation of the drag 
coefficient. The corrected solution on a 5 x 10 3 nodes grid provides the same result as the solution on a 
8 x 10 3 nodes grid. The moderate improvement is due to, the mesh adaptation context that introduces 
a side effect: an adapted mesh with four times the number of vertices of a given adapted mesh is not its 
subdivision by two. Indeed, it has been further improved by the adaptation and provides a better solution 
than a subdivision. We will try to quantify this effect with the next test case. 

Still, the primal and corrected solutions converge to the same value, so that the difference between both 
of them gives and accurate estimation of the remaining numerical error. This second property is particularly 
important in both adaptative and fixed grids context as it informs on how confident we can be in a given 
solution. It also indicates if the grid convergence study should be continued or not. 
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Figure 4. Adapted mesh generated using goal- 
oriented error estimate on the drag for NACA0012 

at Re = 500, M — 0.1 and a = 3°: 8009 vertices. 



Figure 5. Velocity contours around a NACA0012 

at Re = 500, M = 0.1 and a ■ 3°. 


B. 2D Turbulent NACA0012 


The corrector is then validated on a turbulent case with the same geometry, at a Reynolds number of 5 x 10 5 
with the Sapalart-Allmaras turbulence model. The Mach number is fixed to M = 0.7 and the angle of attack 
to a = 1°. Adapted meshes are generated iteratively with f ef lo. a using a hessian-based error estimate based 
on the Mach number field. An example of the resulting meshes and solutions is shown in Figures 6 and 7. 
Additionally, we generate for each mesh its h/2 subdivision and compute the solution on it for comparison. 

Figures 9 and 10 show the convergence of the drag computed with a first and second order scheme for both 
the initial (red) and corrected (black) solution. The blue curve represents for each mesh the drag computed 
on its h/2 mesh subdivision which actual complexity is thus four times larger. Both curves exhibits the same 
trend as the previous laminar example, but it is now possible to quantify how the mesh adaptation improves 
the solution on a mesh with four times more vertices compared to the actual li/2-mesh subdivison. This 
dims the apparent efficiency of the non-linear corrector. In particular we can see that the corrector is very 
efficient with a first order scheme as it yields the same solution as the h/ 2-mesh. It is slightly less efficient 
with a second order scheme but still satisfying. 




Figure 6. Adapted mesh generated using hessian- 

based error estimate on the Mach number field for Figure 7. Velocity contours around a NACA0012 at 
NACA0012 at Re = 5 X 10 5 M = 0.7 and a = 1°: 5394 Re = 5 X 10 5 , M = 0.7 and a = 1°. 
vertices . 


C. 3D inviscid supersonic case: Sonic Boom Prediction Workshop 

The last example is a supersonic case on a realistic aircraft configuration. The considered geometry was 
proposed for the 2 nd AIAA Sonic-Boom Prediction Workshop which aims at minimizing the noise produced 
by supersonic aircrafts on the ground. Different geometries, meshes and configurations were proposed, we 
focus here on the C25D flow through Euler case, see Figures 11 and 12. In order to properly lead the shape 
optimisation process, the noise signature under the aircraft has to be accurate, this is why a grid convergence 
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Figure 8. Convergence of the drag on NACA0012 for initial solution (red) and corrected solution (blue) at 
Reynolds number of 500. 




Figure 9. Convergence of the drag with a first order Figure 10. Convergence of the drag with a second 
scheme on NACA0012 for the initial solution (red), order scheme on NACA0012 for the initial solution 
the corrected solution (black) and the solution on the (red), the corrected solution (black) and the solution 
mesh divided by two (blue) at Reynolds number of on the mesh divided by two (blue) at Reynolds number 

5 x 10 5 . of 5 x 10 5 . 


study is required to determine the numerical error. To this end, mesh adaptation has proven to be very 
efficient 23 and we expect the corrector to provide a second pertinent indication on the convergence of the 
solution. 

Figures 13, 15 and 17 show the pressure signatures (pressure difference with the free stream pressure 
normalized) under the plane computed on tailored grids provided for the workshop. Error bars represent the 
estimated error of each nodal value computed with the nonlinear corrector. We can see that the corrector 
indicates that the error level remains relatively high, even on the 13 million grid, which has been confirmed 
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by the workshop: the grid convergence is not achieved. This is visible on the leading shock which is still 
relatively smooth and the shocks in the aft part of the signature (65 < X < 80) which are still forming. 
Note that the actual subdivision by two of this 13 millions nodes mesh would have yielded a 104 millions 
nodes mesh. 

Figures 14, 16 and 18 shows the same computations using goal-oriented mesh adaptation which have 
been generated to optimally capture the pressure signature. 6 We can se that the error bars are much smaller 
and tend to reduce during the adaptation process. This is especially visible on the mid part of the signature 
(50 < X < 60). It is even more enlightening to compare the error levels predicted by the corrector on the 
tailored grid composed of 13 million nodes in the mid part of the signature to the actual difference between 
the solutions on the tailored grid and adapted grid composed of 5.9 million nodes. The error predicted by 
the corrector on tailored grid matches the difference between the solution on the tailored grid (which is not 
converged) and the solution on the adapted grid (which is converged in this region). 



Figure 11. C25D Geometry of the 2 nd Sonic-Boom 
Prediction Workshop. 



Figure 12. Cut in the volume of the goal-oriented 
adapted mesh and solution on the C25D geometry. 


D. 3D ONERA M6 Wing turbulent Case 

We show here the improvements this method can provide on fixed grids in 3D. We choose to use here the 
ONERA M6 Wing configuration at Mach number M = 0.1, angle of attack a = 3.06° and Reynolds number 
of Re = 11 x 10 6 . A series of four grids were generated using AFLR, 24,25 with a dimensionless wall spacing 
based on the Reynolds number of y + = 3, 2, 1, 0.5. The coarsest mesh is shown in Figure 20. 

We can see in Figure 19 that here, the corrector significantly improves the drag prediction, it yields on 
the 5 x 10 5 nodes mesh the same drag as the solution without correction on the 3 x 10 6 nodes mesh. The 
corrector is even more efficient when the error is not optimally reduced between two grids. We also have 
a good estimation of the quality of our last solution, without the need of any finer mesh: despite the fine 
resolution of the boundary layer, we can be pretty confident on the fact that the solution is not converged 
yet. This is mainly due to the fact that the wake of the wing is poorly discretized, as can be seen in Figure 20. 
This pollutes the computation of the drag. 


IV. Conclusion 

We have seen that the nonlinear correction for inviscid, laminar and RANS flow solutions already shows 
promising results. It improves the overall solution and output functionals, and thus provides a reliable 
estimation of the numerical error. As it does not requires to effectively generate a finer mesh, it is relatively 
cheap to use and can provide a pertinent error estimation even on a single grid. It is particularly efficient 
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on fixed grids and can be used to asses the quality of the mesh refinement that has been done. In a mesh 
adaptation context its apparent efficiency is reduced as the mesh adaptation process efficiently reduces the 
error. Though, it is still a useful tool to asses the convergence of the overall process. However, the transonic 
NACA0012 case showed that the correction is less efficient with a second order scheme than with a first 
order scheme, suggesting that the MUSCL extrapolation requires an additional particular treatment. 
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Tailored-Adaptation Error Bars on C25D Geometry 3.4M - R/L=l - Phi=0 



X 

Figure 13. Pressure difference under the plane 
for tailored grid of 3.4 million vertices (red line) 
with estimated error provided by nonlinear correc¬ 
tor (black error bars). 


Tailored-Adaptation Error Bars on C25D Geometry 6.3M - R/L=l - Phi=0 
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Goal-Oriented Adaptation Error Bars on C25D Geometry 1.5M - R/L=l - Phi=0 



X 

Figure 14. Pressure difference under the plane for 
goal-oriented adapted grid of 1.5 million vertices 
(red line) with estimated error provided by nonlin¬ 
ear corrector (black error bars). 


Goal-Oriented Adaptation Error Bars on C25D Geometry 3.1M - R/L=l - Phi=0 
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Figure 15. Pressure difference under the plane 
for tailored grid of 6.3 million vertices (red line) 
with estimated error provided by nonlinear correc¬ 
tor (black error bars). 


Figure 16. Pressure difference under the plane for 
goal-oriented adapted grid of 3.1 million vertices 
(red line) with estimated error provided by nonlin¬ 
ear corrector (black error bars). 


Tailored-Adaptation Error Bars on C25D Geometry 13.0M - R/L=l - Phi=0 



Goal-Oriented Adaptation Error Bars on C25D Geometry 5.9M - R/L=l - Phi=0 



X 
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Figure IT. Pressure difference under the plane 
for tailored grid of 13 million vertices (red line) 
with estimated error provided by nonlinear correc¬ 
tor (black error bars). 


Figure 18. Pressure difference under the plane for 
goal-oriented adapted grid of 5.9 million vertices 
(red line) with estimated error provided by nonlin¬ 
ear corrector (black error bars). 
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Figure 19. Convergence of the drag on ONERA 
M6 Wing for initial solution (red) and corrected 
solution (green) at Reynolds number of 11 x 10 6 
and angle of attack a = 3.06°. 


Figure 20. Coarsest M6 Wing mesh generated with 
a fixed structured boundary layer and frontal ap¬ 
proach with AFLR: 110009 vertices. 
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